# TODO: Add comment
# 
# Author: broken
###############################################################################
library(raster)
library(base)
library(gstat)
library(RPostgreSQL)
dr=dbDriver("PostgreSQL")

path= "C:/Users/forever/Desktop/MYD01/"

con=dbConnect(dr,dbname="fimo",host="192.168.0.4",port=5432,user="rasdaman",password="rasdaman")
timeQuery="SELECT aqstime FROM apom.satresampmyd04 where filename like '%"
selectStr="SELECT stationid,name,st_asText(location),avg(pm25) as avg_pm FROM apom.grdpmhiscem_data,public.grdstation WHERE grdpmhiscem_data.stationid = grdstation.id AND grdstation.id IN (1,2,4,5) AND  aqstime between '"
andStr="' and '"
groupStr="'  group by stationid,name,location having avg(pm25) > 0"

year_array=c(2010,2011,2012,2013,2014)
for (k in 1:5) {
	inputPath = paste(path,year_array[k],sep="")
	listFile = list.files(path = inputPath,pattern="10km _og.tif$",full.names = FALSE,recursive=TRUE)
	listFile=substr(listFile,1,40)
	listSize=length(listFile)
	
	listTime=c(1:listSize)
	
	## sCor=sRMSE=sRE=c(1:listSize)
	## 
	## oCor=oRMSE=oRE=c(1:listSize)
	
	uCor=uRMSE=uRE=c(1:listSize)
	totalRow=c(1:listSize)
	
	for(i in 1:listSize){
		fileName=listFile[i]
		
		ogFile=paste(inputPath,"/",fileName,".hdf_CB_10km _og.tif",sep="")
		## skFile=paste(inputPath,"/",fileName,"_sk.tif",sep="")
		## okFile=paste(inputPath,"/",fileName,"_ok.tif",sep="")
		ukFile=paste(inputPath,"/",fileName,"_uk.tif",sep="")
	
		ogRaster=raster(ogFile)
		## skRaster=raster(skFile)
		## okRaster=raster(okFile)
		ukRaster=raster(ukFile)
	
		qr=paste(timeQuery,fileName,"%'",sep="")
		resultSet=dbSendQuery(con,qr)
		record=fetch(resultSet,n=1)
		asqTime=substr(record[1,1],1,30)
		myTime=strptime(asqTime,format="%Y-%m-%d %H:%M:%S")
	
		#add 7h
		myTime=myTime+25200
		
		startTime=myTime-1800
		endTime=myTime+1800
		
		
		selectQuery=paste(selectStr,startTime,andStr,endTime,groupStr,sep="")
		pmRs=dbSendQuery(con,selectQuery)
		pmRecord=fetch(pmRs,n=-1)
		numberRow=nrow(pmRecord)
		if(numberRow>0){
			ogValueArray=c(1:numberRow)
			skValueArray=c(1:numberRow)
			okValueArray=c(1:numberRow)
			ukValueArray=c(1:numberRow)
		
			for(j in 1:numberRow){
				locationStr=pmRecord[j,3]
				index=regexpr(" ",locationStr)[1]
				xLocation=substr(locationStr,7,index-1)
				yLocation=substr(locationStr,index+1,nchar(locationStr)-1)
				
				ogValue=extract(ogRaster,cbind(as.numeric(xLocation),as.numeric(yLocation)))
				## skValue=extract(skRaster,cbind(as.numeric(xLocation),as.numeric(yLocation)),buffer=25000,fun=mean)
				## okValue=extract(okRaster,cbind(as.numeric(xLocation),as.numeric(yLocation)),buffer=25000,fun=mean)
				ukValue=extract(ukRaster,cbind(as.numeric(xLocation),as.numeric(yLocation)),buffer=25000,fun=mean)
				
				ogValueArray[j]=ogValue
				## skValueArray[j]=skValue
				## okValueArray[j]=okValue
				ukValueArray[j]=ukValue
			
			}
			pmRecord$ogColumn=ogValueArray
			## pmRecord$skColumn=skValueArray
			## pmRecord$okColumn=okValueArray
			pmRecord$ukColumn=ukValueArray
			
			## s_cor=cor(pmRecord$avg_pm,pmRecord$skColumn)
			## s_rmse=sqrt(sum((pmRecord$skColumn-pmRecord$avg_pm)^2)/numberRow)
			## s_re=sum(abs(pmRecord$skColumn-pmRecord$avg_pm)/pmRecord$avg_pm)/numberRow
			## 
			## sCor[i]=s_cor*s_cor
			## sRMSE[i]=s_rmse
			## sRE[i]=s_re*100
			totalRow[i]=numberRow
			
			## 
			## o_cor=cor(pmRecord$avg_pm,pmRecord$okColumn)
			## o_rmse=sqrt(sum((pmRecord$okColumn-pmRecord$avg_pm)^2)/numberRow)
			## o_re=sum(abs(pmRecord$okColumn-pmRecord$avg_pm)/pmRecord$avg_pm)/numberRow
			## 
			## oCor[i]=o_cor*o_cor
			## oRMSE[i]=o_rmse
			## oRE[i]=o_re*100
			
			u_cor=cor(pmRecord$avg_pm,pmRecord$ukColumn)
			u_rmse=sqrt(sum((pmRecord$ukColumn-pmRecord$avg_pm)^2)/numberRow)
			u_re=sum(abs(pmRecord$ukColumn-pmRecord$avg_pm)/pmRecord$avg_pm)/numberRow
			
			uCor[i]=u_cor*u_cor
			uRMSE[i]=u_rmse
			uRE[i]=u_re*100
			
			outFile=paste(inputPath,"/",fileName,"_validate.csv",sep="")
			write.csv(pmRecord,file=outFile)
			print(paste("Complete ",fileName,sep=" "))
		
		}else{
			print(paste(fileName,"have no data ",sep=" "))
		}
		
	}
	result_dataframe=data.frame(listFile,totalRow,uCor,uRMSE,uRE)
	outFile=paste(inputPath,"/","mod_year_4np.csv",sep="")
	write.csv(result_dataframe,file=outFile)
	print(paste("Finish",year_array[k]))
	
}

print("Finish All")




